Rozważamy dwie masy punktowe, M i m, orbitujące po okręgach dookoła wspólnego środka masy. Taki uproszczony problem 2 ciał (ogólnie ruch jest po elipsach, paraboli lub hiperboli) stanowi wyjście do dyskusji sił pływowych i najprostszego rezonansu: orbitalno-rotacyjnego. Najbliższym przykładem jest system Ziemia-Księżyc.

Solve[{M Ω^2 dM==G M m/d^2,m Ω^2 dm==G M m/d^2,dM+dm==d},{Ω,dM,dm}] (* Z warunki równości przyspieszenia grawitacyjnego i odśrodkowego *)

rozw=%//Last  (* Dwa rozwiązania na częstość keplerowską Ω odpowiadają rotacji przeciwnie i zgodnie do ruchu wskazówek zegara. Bierzemy to ostatnie. *)

Ω^2/.rozw (* Częstość keplerowska zależy od SUMY mas *)

{dM,dm,dM+dm}/.rozw/.m->q M//Simplify (* Można wprowadzić stosunek mas q=m/M. Użycie różnej wielkości liter podpowiada, że często M>>m, ale wszystkie wzory muszą wyjść symetryczne zw względu na zamianę m  <> M. Czyli na zamianę q->1/q *)

%/.q->1/q

Simplify[%]

Inny sposób wyprowadzenia tego samego, wprost z równań dynamiki Newtona.

Clear[m,M,T,R,d]

(* Równanie ruchu Newtona *)

eq1a=m x1''[t]==
-G m M (x1[t]-x2[t])/Sqrt[(x1[t]-x2[t])^2+(y1[t]-y2[t])^2]^3

eq1b=m y1''[t]==
-G m M (y1[t]-y2[t])/Sqrt[(x1[t]-x2[t])^2+(y1[t]-y2[t])^2]^3

eq2a=M x2''[t]==
-G m M (x2[t]-x1[t])/Sqrt[(x2[t]-x1[t])^2+(y2[t]-y1[t])^2]^3

eq2b=M y2''[t]==
-G m M (y2[t]-y1[t])/Sqrt[(x2[t]-x1[t])^2+(y2[t]-y1[t])^2]^3

sys={eq1a,eq1b,eq2a,eq2b};

(* Podstawienie do równań ruchu po okręgu *)

sys/.{
x1->(M/(m+M) d*Cos[Ω*#]&),
y1->(M/(m+M) d*Sin[Ω*#]&),
x2->(-m/(m+M)d*Cos[Ω*#]&),
y2->(-m/(m+M)d*Sin[Ω*#]&)
};

(* Sprawdzenie wyniku *)

Simplify[%,d>0]

%/.{Ω->Sqrt[ G (m+M)/d^3]}

m=1

M=3

G=d=1

Ω=Sqrt[G (M+m)/d^3]

Tmax=0.85 2Pi/Ω

SetOptions[ParametricPlot,PlotRange->{{-1,1},{-1,1}},Ticks->{{{-1,"-d"},{-M/(m+M),"-M/(m+M)"},{M/(m+M),"M/(m+M)"},{1,"d"}},{{-1,"-d"},{-m/(m+M),"-m/(m+M)"},{m/(m+M),"m/(m+M)"},{1,"d"}}}];

w=Block[{T=Tmax},ParametricPlot[
{M/(m+M){ Cos[Ω t],Sin[Ω t]},-m/(m+M){Cos[Ω t],Sin[Ω t]}},{t,0,T},
Epilog->{Disk[M/(m+M){ Cos[Ω T],Sin[Ω T]},0.02],Disk[-m/(m+M){ Cos[Ω T],Sin[Ω T]},0.05]},PlotStyle->{{Red,Thick},{Blue,Thick}}]]

wM=Block[{T=Tmax},ParametricPlot[{M/(m+M){ Cos[Ω t],Sin[Ω t]}+m/(m+M){Cos[Ω t],Sin[Ω t]}},{t,0,T},
Epilog->{Disk[{ Cos[Ω T],Sin[Ω T]},0.02],Disk[{0,0},0.05]},PlotStyle->{{Red,Thick},{Blue,Thick}}]]

wm=Block[{T=Tmax},ParametricPlot[-{M/(m+M){ Cos[Ω t],Sin[Ω t]}+m/(m+M){Cos[Ω t],Sin[Ω t]}},{t,0,T},Epilog->{Disk[{ 0,0},0.02],Disk[-{ Cos[Ω T],Sin[Ω T]},0.05]},PlotStyle->{{Blue,Thick}}]]

GraphicsGrid[{{w,wM,wm}}]

Clear[T]

Manipulate[
GraphicsGrid[
{
{
ParametricPlot[
{M/(m+M){ Cos[Ω t],Sin[Ω t]},-m/(m+M){Cos[Ω t],Sin[Ω t]}},{t,0,T},
Epilog->{Disk[M/(m+M){ Cos[Ω T],Sin[Ω T]},0.02],Disk[-m/(m+M){ Cos[Ω T],Sin[Ω T]},0.05]},PlotStyle->{{Red,Thick},{Blue,Thick}}],
ParametricPlot[{M/(m+M){ Cos[Ω t],Sin[Ω t]}+m/(m+M){Cos[Ω t],Sin[Ω t]}},{t,0,T},
Epilog->{Disk[{ Cos[Ω T],Sin[Ω T]},0.02],Disk[{0,0},0.05]},PlotStyle->{{Red,Thick},{Blue,Thick}}],
ParametricPlot[-{M/(m+M){ Cos[Ω t],Sin[Ω t]}+m/(m+M){Cos[Ω t],Sin[Ω t]}},{t,0,T},Epilog->{Disk[{ 0,0},0.02],Disk[-{ Cos[Ω T],Sin[Ω T]},0.05]},PlotStyle->{{Blue,Thick}}]
}
}],{{T,Tmax},-2Pi/Ω,2Pi/Ω}]